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Binary black-hole coalescence is treated semi-analytically by a novel approach. Our prescription 
employs the conservative Skeleton Hamiltonian that describes orbiting Brill-Lindquist wormholes 
(termed punctures in Numerical Relativity) within a waveless truncation to the Einstein field equa- 
tions [G. Faye, P. Jaranowski and G. Schafer, Phys. Rev. D 69, 124029 (2004)]. We incorporate, 
in a transparent Hamiltonian way and in Burke-Thorne gauge structure, the effects of gravita- 
tional radiation reaction into the above Skeleton dynamics with the help of 3.5PN accurate angular 
momentum flux for compact binaries in quasi-circular orbits to obtain a Semi-Analytic Puncture 
Evolution to model merging black-hole binaries. With the help of the Taylor T4 approximant at 
3.5PN order, we perform a first-order comparison between gravitational wave phase evolutions in 
Numerical Relativity and our approach for equal-mass binary black holes. This comparison reveals 
that a modified Skeletonian reactive dynamics that employs flexible parameters will be required to 
prevent the dephasing between our scheme and Numerical Relativity, similar to what is pursued in 
the Effective One Body approach. A rough estimate for the gravitational waveform associated with 
the binary black-hole coalescence in our approach is also provided. 

PACS numbers: 04.30.Db, 04.25.Nx 



I. INTRODUCTION 

Binary black hole mergers are the most plausible 
sources of gravitational radiation for the currently op- 
erational and planned ground-based laser interferometric 
gravitational-wave (GW) detectors like GEO-LIGO and 
VIRGO. Further, coalescing binary black holes (BBH) 
are considered to be the workhorse sources to realize GW 
astronomy with the proposed space-based GW interfer- 
ometer LISA. Recent advances in Numerical Relativity, 
achieved by several independent groups, should enable 
its practitioners to create accurate and 'exact' numer- 
ical relativity based GW search templates for the late 
stages of binary black hole coalescence [lj. However, it 
is also desirable to develop semi-analytic prescriptions, 
based on well defined assumptions, capable of describing 
the most relativistic part of the binary black hole coales- 
cence, namely the last few cycles of the inspiral, plunge 
and the subsequent merger. The two main reasons to de- 
velop such descriptions are as follows. It is reasonable to 
imagine that semi-analytic approaches, based on certain 
assumptions, should be able to probe the purely numer- 
ical approaches in order to extract the 'physics' behind 
binary black hole coalescence. The second reason is re- 
lated to the GW data analysis requirement of creating 
huge banks of accurate GW templates to cover, rather 
densely, the binary black hole parameter space. 

In the literature and to best of our knowledge, there ex- 
ists only two approaches that can semi-analytically model 
the late stages of non-spinning compact binary coales- 
cence. Historically, the first prescription is the so-called 
hybrid approach, developed by Will and his collabora- 
tors, which re-sums exactly the 'Schwarzschild terms' in 
the post-Newtonian(PN) accurate equations of motion 
for compact binaries and treats the symmetric mass ra- 
tio terms as additional corrections Hi. The second and 



widely employed approach is the so-called Effective One- 
Body (EOB) scheme, developed by Damour and collab- 
orators [1, The EOB scheme requires available PN 
accurate results for the conservative orbital dynamics 
and GW luminosity associated with comparable mass 
compact binaries moving along quasi-circular orbits ex- 
tractable from Refs. [5|. In the EOB approach, one pack- 
ages these PN accurate results in a certain re-summed 
form with the hope of extending the validity of PN ac- 
curate results to strong- field (and fast-motion) regimes. 
The EOB package allows PN accurate dynamics of a com- 
parable mass compact binary, having masses mi and m,2, 
to be mapped onto an essentially geodesic dynamics of 
a test-particle of mass /i = m\ raij(ra\ + 1712) moving in 
an 77 = mi 7712/ (mi + m-i) 2 deformed Schwarzschild ge- 
ometry. Further, reactive dynamics is included with the 
help of certain Pade re-summed radiation reaction force. 
It is important to note that in the hybrid and EOB ap- 
proaches, describing late stages of binary black hole coa- 
lescence, essentially contain a single black hole and some 
undesirable features of these schemes were pointed out in 
Refs. @ and 0. 

In this paper, we present a prescription that can semi- 
analytically model the late stages of coalescing non- 
spinning binary black holes (BBH). Our approach em- 
ploys the 'Skeleton' solution to the Einstein field equa- 
tions, detailed in Ref. Q, to describe the BBH conserva- 
tive dynamics. We augment the above mentioned BBH 
conservative dynamics by incorporating radiation reac- 
tion effects in a fully Hamiltonian way with the help of 
the 3.5PN accurate GW luminosity, available in Refs. [f|. 
This is how, in a nutshell, we provide a semi-analytical 
description for the late stages of BBH coalescence. The 
fact that the Skeleton Hamiltonian is constructed in the 
Arnowitt-Deser-Misner (ADM) 3+1 approach to General 
Relativity |9| makes our scheme very close to the BBH dy- 
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namics, analyzed by several Numerical Relativity groups. 
The Skeleton Hamiltonian contains, in its static part, the 
potential of the well-known Brill-Lindquist BBH initial 
value configuration [l(|[D]]- Further, it allows black holes 
to orbit each other along conservative orbital configura- 
tions. We emphasize that the above mentioned conser- 
vative orbital evolution obey exactly only a truncated 
version of Einstein's field equations Q. In the Skele- 
ton approach for black holes in circular orbits, Ref. [1] 
also demonstrated an intriguing property; the ADM mass 
(identified through the 3-metric) coincides with the 
Komar mass (identified through the pure time compo- 
nent of the 4-metric) (l2| . We note that the equality of 
ADM and Komar masses is a general property of station- 
ary solutions to the Einstein field equations for isolated 
systems. This is, in our view, an appealing feature of the 
Skeleton scheme. 

In Numerical Relativity simulations that employ punc- 
tures to model black holes, it is customary to use Brill- 
Lindquist configurations as initial data [1 31 ] - Therefore, 
the unique construction of the Skeleton Hamiltonian for 
BBH, based on certain well defined assumptions, and 
its present extension should allow us to probe the var- 
ious aspects of Numerical Relativity based BBH simula- 
tions involving punctures. This prompted us to name our 
approach as Semi- Analytic Puncture Evolution (SAPE). 
However, it should be noted that the recent Numerical 
Relativity simulations involving black holes do not im- 
pose conformal flat condition in their evolutions, required 
by our SAPE. 

The organization of this paper is as follows. In the 
next section, we briefly summarize the Skeleton approxi- 
mation and the conservative binary black hole dynamics 
it defines, detailed in Ref. Section [II] also includes 
a fully consistent Hamiltonian way of incorporating the 
effects of radiation reaction on to the conservative Skele- 
ton dynamics. The resulting binary black hole dynamics 
and its implications are explored in Section IIIII An es- 
timate for the complete gravitational wave polarizations 
associated with the binary black hole coalescence in re- 
active Skeleton dynamics, in the 'restricted waveform' 
approach, is presented in Section IIV1 Our conclusions 
and future directions are available in Section fVl 



A. Conservative Skeletonian Binary Black Hole 
Dynamics 

In the ADM formulation of General Relativity, the 
spacetime line element in the (3 + 1) decomposed form is 
given by 

ds 2 = -a 2 c 2 dt 2 + ^{dx l + I3 l cdt)(dx ] + (3 1 cdt) 

(ij = 1,2,3), (I) 

where a is the lapse function, j3 l the shift vector, 7y the 
induced metric on a three-dimensional spatial slice S(i), 
parametrized by the time coordinate t, and c is the veloc- 
ity of light. The canonical conjugate to the 3-metric 7^ 
is given by 7r y c 3 /167rG, where G is the Newtonian gravi- 
tational constant and it is a tensor density of weight +1. 
Further, 7^ and tt 13 satisfy the Hamiltonian and momen- 
tum constraints [9]. On any three-dimensional spatial 
slice £(i), the above four constraints may be written as 

^-^('W-5' 2 )^?("'" c2 

\l/2 

(2a) 

-29 i7 r| + K kl d akl = ^ ]T p a iS a , (2b) 

a 

where R, 7, 7 y stand respectively for the curvature 
scalar, the determinant, the inverse metric associated 
with the 3-metric 7^ . The linear momentum of point- 
mass a (a=l,2) is denoted by p ai and its bare mass by m a . 
The partial derivative with respect to the hypersurface 
coordinates is denoted by di — d/dx 1 . The Dirac delta 
function S a — S(x l — x l a (t)) is defined by 5 (x l —x l a (tj) = 
if x 1 ^ x l a {t) and / d 3 x d(x l — x l a (t)) = 1, where x l a {t) is 
the space coordinate of point-mass a at time t. 

The ADM coordinate conditions that generalize the 
isotropic Schwarzschild metric read 



7« = ^1 + g< 

7T = U , 



■ TT 



(3a) 
(3b) 



II. BINARY PUNCTURE DYNAMICS IN 
SKELETON APPROXIMATION 



where is a function vanishing at spatial infinity. The 
tranverse-traceless (TT) part of the metric 7^ with re- 
spect to the Euclidean 3-metric is denoted by hf^ 
and the Einstein summation convention is also employed 
in the above equation, 7r M = 0. Taking into account the 
above gauge condition for ir lJ , the following decomposi- 
tion can be achieved: 



In this section, we briefly summarize the construction 
of the conservative binary black hole Skeleton Hamilto- 
nian, detailed in Ref. [8(. Afterwards, we explain our 
prescription for including the effects of GW radiation re- 
action into the Skeletonian conservative dynamics. 



(4) 



where n lJ denotes the longitudinal part of tt 11 which may 
be expressed as 



fr« = dnr j + d 3 ix l - -6 t] d k ir k , 



(5) 
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where it 1 is a vector in the flat 3-space. Further, the TT 
part of 7r u , or rather 7r^ T c 3 /l 6irG, provides the canonical 
conjugate to hj? . 

The Hamilton functional for non-spinning binary black 
hole system is given by 



H 



16ttG 

,TT ij 
l ij j " TT 



d 3 x Acj) 



x a ! Pai i 



(6) 



where A is the Laplacian in the 3-dimensional flat space. 
While dealing with the dynamics of binary black-hole 
systems, it is highly advantageous to introduce a Routh 
functional of the form, 



R [x l a , Pai , hj?, d t hj?} = H - ^ J 



d 3 x Tr^ T d t hJ T . 



(7) 



and as expected it describes the total (initial) energy be- 
tween two uncharged Brill-Lindquist black holes. 

The Skeleton approach to General Relativity for deal- 
ing with BBH space-times, detailed Ref. [8J, requires both 
the conformal flat condition/truncation for the spatial 3- 
metric and hereof together with Eq. (|3bj) the maximal 
slicing condition. These prescriptions, valid for all times, 
are given by 



1 



Hi = C 1 + §<W% > 
7r y 7y =2^ ij Kij =0. 



(11a) 
(lib) 



is the extrinsic curvature of £(t) 



In above equations, Ki 3 

and is related to tt^ by7r« = -7 1 / 2 (7 i V m -7 i V m )^m- 
Under the above conditions, the momentum constraint 
equations for the Einstein theory become 



8ttG 



E 



Pai5 a ■ 



(12) 



This functional is a Hamiltonian for the particle (point- 
mass) degrees of freedom, and a Lagrangian for the in- 
dependent gravitational field degrees of freedom. The 
truncation implemented in Ref. [8j, i.e. to put hj^ = 0, 
is called the conformal flat condition and it results in 
H = R{xi,p ai ). 

It is possible to show that the introduced source model, 
namely S a terms in Eqs. @, can produce certain bi- 
nary black-hole spacetime. Indeed, it was demonstrated 
in Refs. [ll| that one obtains the Brill-Lindquist initial 
value solution for uncharged black holes at rest, available 
in Ref. fioll. Therefore, let us first summarize Refs. [TT[ | . 
In Refs. [ll| , one imposes hJ } T — and p ai — (and this 
leads to it 

nian constraint becomes 



0). Under these restrictions, the Hamilto- 



Acj> 



16ttG 



m a S a . 



(8) 



The exact (and unique) solution to the above equation, 
detailed in Ref. [ll| and obtained with the help of dimen- 
sional regularization, reads 



4G ^ a a 
— 1^ T ' where 



c 
m 



(9a) 



a a 



rrib 



2c 2 r ab /G 



c 2 r ab 
G 

2-i 1/2 



c 2 r ab /G 



(9b) 



and r a denotes the Euclidean distance between x 1 and x l a 
and r ab the Euclidean distance between x l a and x\. 

The Brill-Lindquist potential function, computed in 
Refs. [11| . reads 



Hbl = ("i + 012) c 2 = (mi + m 2 ) c 2 - G 



ri2 



(10) 



The solution of this equation is constructed under the 
condition that irf is of purely longitudinal form 



t} = by, + aw - -SijW 



(13) 



where Vi is a vector in the flat 3-space. The difference to 
the way 7r y is defined, via Eq. (JSJ), should be noted. In 
the Skeleton approach, nj consists only of its longitudinal 
part and is defined in terms of Vi via Eq. (TT3|) . The 
Hamiltonian constraint in the Skeleton scheme becomes 



16ttG v ^ 
~~ ^ ^ (1 



Pa 



(1 



c° 

1/2 



m a 8g 



1+ 



(14) 



Further, one implements a truncation to the numerator 
of the first term, namely 7r^7rj, in the following manner 



16ttG 



J^PajVjSa. 



(15) 



This is the other crucial truncation in addition to the 
conformal flat truncation required in the Skeleton ap- 
proach. This truncation introduces errors of the same 
order as the conformal flat truncation. Therefore, in a 
post-Newtonian setting, the errors are of the 2PN order. 
With the aid of the ansatz 



n 2 /L^i 



(16) 



which is inspired by the delta-function sources in 
Eq. (fl"4"|) , after the insertion of Eq. (fT5"|) . the resulting 
energy constraint equation can be solved in the sense of 
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dimensional regularization, detailed in Ref. @. This re- 
sults in two coupled algebraic equations given by 



1 + A S 



1 P 2 a/(m 2 a C 2 ) 
PaiVai/c 



(17) 



where a = 1,2, 6 ^ a and A stands for G/(2c 2 ). The 
a a -parameters generalize the rest energies of interacting 
black holes at rest, given in Eq. (|9bj). to orbiting black 
holes. 

With these inputs, the Skeleton Hamiltonian for BBH 
becomes 

4 r 

H ^ = -^-Gj ^A0 = c 2 ^a a . (18) 
The Hamilton equations of motion, as usual, read 

X =™ v=~— (19) 
a -dp a > Pa dx a [U) 

For a binary system in the center-of-mass frame, we have 

Pi = P2 = P , (20a) 
r = X! - x 2 , (20b) 
r 2 = (xi - x 2 ) ■ (x-i - x 2 ) . (20c) 

Further, we will employ the following convenient di- 
mensionless quantities 



t 



J 



tc A 
Gm : 
Jc 



r c 
Gm 
P 



P H Sk 

7 P = ! HSk = n 

flC [ICS 



, Pr = ~, p 2 =pl+?/f\ 
fl C 



(21a) 
(21b) 



Gm n 

where J = r x p and p r = p ■ r/r. 

For the purpose of present investigation and using 
Eqs. (5.1), (5. 2) and (5.4) in Ref. Q, we write the Skele- 
ton Hamiltonian 7isk in the following way 



n 



Sk 



i>i = 1 



2 f [ ipi + ip 2 — 2j , where 

^V 2 (p 2 r + w) 



(22a) 



X- 



4f 1P2 
(8p 2 + 7l 



1 



if 



8f 2 V>2 7 



lp 2 = l + 



X- 



4r? 2 {£ + £ 

x\¥i 



1/2 



(22b) 



4f ipi 
(8p 2 + 7| 



1/2 



8r 2 V-i 7 



(22c) 



and the symbols X-/+ are defined to be \- = 
(1 - ^/T - 477) and x+ = (1 + 7^477). 

The conservative Skeleton Hamiltonian has the follow- 
ing features and properties. The Skeleton Hamiltonian is 
exact in the test-body limit, where it describes the mo- 
tion of a test particle in the Schwarzschild spacetime. It 
is also identical to the 1PN accurate Hamiltonian describ- 
ing the compact binary dynamics in General Relativity. 
Further, as explained earlier, when the point particles 
are at rest, the Brill-Lindquist initial value solution is 
reproduced. It is remarkable that the Skeleton Hamil- 
tonian allows a post-Newtonian expansion in powers of 
1/c 2 to arbitrary orders. This is remarkable because of 
the observation that the 3PN approximate binary black 
hole Hamiltonian does not exist, if one only imposes the 
conformal flat condition without invoking Eq. (fT5|) 
Therefore, the Skeleton Hamiltonian describes the evo- 
lution of punctures under both conformal flat conditions 
for the 3-metric and analyticity conditions for the Hamil- 
tonian via PN expansions (and, of course, gravitational 
radiation reaction is not incorporated) . The Skeleton dy- 
namics is not isolated from the metric field rather it is 
embedded in the field [recall that the Skeleton Hamilto- 
nian originates from 7^]. Furthermore, with the aid of 
similar truncations, it is also possible to calculate a and 
P l [8j and, as expected, the demonstration of the equal- 
ity of ADM and Komar masses for BBH in circular orbits 
in the Skelton approach requires explicit expressions for 
a and f3 l . These arguments imply that in the Skelton 
scheme for BBH, a global spacetime picture exists. 



B. Tackling the effects of radiation reaction 

This subsection details how we incorporate radiation 
reaction effects into the above described conservative bi- 
nary black hole dynamics. Our prescription is heavily 
influenced by the way the reactive 2.5PN Hamiltonian 
affects the conservative 2PN accurate compact binary 
dynamics. Therefore, let us first take a closer look at 
the way reactive contributions to TL change Hamiltonian 
equations of motion in the post-Newtonian approxima- 
tion with the help of 2.5PN accurate Tt, available in 
Refs. [HI, EH]. It is convenient (and usual) to split the 
fully 2.5PN accurate H as 



H- 



2.5PN 



H: 



2PN 



H. 



2.5PN ! 



(23) 



where -ff^PN defines the 2PN accurate conservative or- 
bital dynamics [15]. The explicitly time dependent re- 
active Hamiltonian, whose first contribution appears at 
2.5PN order, is denoted by H, 



2.5PN' 



In terms of the un- 
sealed position and momentum vector components in the 
center-of-mass frame, and following Ref. [14], it reads 



_ 2G d 3 Q 



(24) 
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where Qij = [i {vifj — r 2 <%/3) stands for the Newtonian 
accurate mass-quadrupole of a compact binary. How- 
ever, for the purpose of our investigation, it is convenient 
to consider the following physically equivalent reactive 
Hamiltonian at the 2.5PN order: 



H, 



2.5PN 



Gfi d b Qij 
5c? dt 5 



(25) 



The above two reactive Hamiltonians differ by a total 
time derivative and thus belong to two different coordi- 
nate conditions [it is only for simplicity that dynamical 
variables are still denoted by same symbols in Eqs. 
and (|25p]. The reactive Hamiltonian given by Eq. 
is valid under the ADM coordinate conditions, whereas 
Eq. (f2"5|) is for its counterpart in the Burke-Thorne gauge 
PJ5| (notice that in the Burke-Thorne gauge, the reactive 

■ TT 



part of hlj is zero as explained in Ref. [17(). However, 
a change of coordinates at the reactive order, in compar- 
ison with the conservative part, does not influence any 
observable quantities. 

With the help of Eq. (|S5]), we define H2.5PN = 
H\ 5PN //J c 2 , and obtain, rather symbolically for the 2PN 
conservative part, the following fully 2.5PN accurate 
Hamiltonian equations of motion in polar coordinates 
and in dimcnsionlcss variables as 



dr 
di 
# 
dt 

dp r 

di 
d± 
di 



dU 



2PN 



dp r 

dj 

^^2PN 



dr 

d^2.5PN _ 



fPr 



dd 



(26a) 
(26b) 
(26c) 
(26d) 



where Ti^PN = ■^2Pn/' uc anc ^ the explicit expressions 
j are gr 

16 77 p r 



for fp r and f~- are given by 



15 f 5 
5 f 5 



4f + 3p 2 r f 2 + 18 j 



9p 2 r f 2 + 6f -2f 



(27a) 
(27b) 



As mentioned earlier, in this investigation we focus on 
black hole binaries inspiralling along quasi-circular orbits 
and we prescribe the quasi-circularity by imposing p r = 
on the reactive contributions to Hamiltonian equations of 
motion. This implies that Eqs. (|27|) becomes 



0. 



32 77 

■5" ff/2 



(28a) 
(28b) 



A close look at the above equation for /- that defines 

dj /dt reveals that it is identical to the dominant contri- 
bution to the far-zone angular momentum flux associated 



with comparable mass binaries inspiralling along quasi- 
circular orbits. This observation helped us to provide 
the following prescription to incorporate radiation reac- 
tion on the conservative Skeleton Hamiltonian, given by 
Eqs. (EU). 

We recall that the Skeleton Hamiltonian providing the 
conservative part of the binary black hole dynamics is 
PN independent and therefore can be employed to treat 
close BBH configurations, provided we accept the valid- 
ity of underlying assumptions. Therefore, in order to 
introduce gravitational radiation reaction effects, we in- 
voke the most PN accurate expression for far-zone an- 
gular momentum flux associated with comparable mass 
compact binaries inspiralling along quasi-circular orbits 
to provide an evolution equation for j. This is how we 
incorporate radiation reaction into Eqs. (|22[) and the dy- 
namics it defines. The fully 3.5PN accurate expression 
for far-zone angular momentum flux for compact binaries 
in quasi-circular orbits, denoted by Lj and extractable 
from Ref. [f| , reads 



32 



44711 



1247 35 



12 



336 
9271 



65 



9072 



-5bT' 7+ 18 7? 



Attv 3 
"8191 



672 



583 



24 

1712 



nv 5 + 



6643739519 16 



105 
775 
324 
214745 

1728 



7 + 



69854400 
134543 41 



7776 



7Z n V 



94403 



1712 



In (4u) 



105 
193385 



3024 
16285 



504 



3024 



' if 



(29) 



where v = lu 1 / 3 and uj = d<p/di is the conservative dimen- 
sionless orbital frequency. The above equation is derived 
by noting that for circular orbits, the far-zone angular 
momentum and energy fluxes are related by Le = ^Lj. 

Therefore, the Skeleton dynamics representing dynam- 
ics of comparable-mass non-spinning coalescing black 
hole binaries, in the scaled polar coordinates, can be sum- 
marized as 



(30a) 
(30b) 
(30c) 
(30d) 



dr 




di ~ 


dp r 


dej) 


cWsk _ 


di ~ 


dj 


dp r 


cWsk 


di 


df 


dl = 




di 





where Lj, the 3.5PN accurate far-zone angular momen- 
tum flux, is given by Eq. (|29p . It is important to note 
that, in our prescription, Cj appearing in the PN accurate 
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expression for Lj via v = ili 1 ^ 3 is provided by the Skele- 
ton dynamics [in other words, by Eq. (|30b[) ]. The com- 
putational details required to compute right hand sides 
of Eqs. (f3T)|) are provided in the Appendix lAl 

We have already argued that the Skeleton Hamilto- 
nian, given by Eq. (|22|) . does provide a certain conser- 
vative dynamics for two punctures orbiting each other. 
Therefore, Eqs. (|30p represent a prescription to im- 
plement semi- analytically puncture evolution to model 
merging binary black holes with reasonable precision. It 
should be noted that the use of Eq. ([29| implies that the 
adiabatic inspiral occurs along quasi-circular orbits. 

We observe that the arguments that allowed us to em- 
ploy PN accurate far-zone angular momentum flux to 
incorporate the effects of radiation reaction into our con- 
servative (and PN independent) Hamiltonian is fairly el- 
egant compared to somewhat heuristic arguments em- 
ployed by the advocates of the EOB scheme (see Sec. Ill 
in Ref. [3). Further, the EOB approach usually re- 
quires re-summed Pade estimates for dj/dt that also re- 
quire PN accurate far-zone angular momentum flux, to 
include effects of radiation reaction on the conservative 
EOB Hamiltonian. We have experimented with various 
Pade estimates for dj/dt in our Skeletonian approach. 
However, the final value of the Skeletonian j, when u) 
reaches its maximum value, turned out to be rather in- 
sensitive to the choice of Pade or Taylor approximant 
for dj/dt. Therefore, in this paper, we employed only 
the usual 3.5PN (Taylor) accurate expression for far-zone 
angular momentum flux, in terms of v = a) 1 / 3 , given by 
Eq. (HU). 

In the next section, we describe how we numerically 
implement the reactive Skeleton dynamics and explore 
its various facets. 



III. NUMERICAL EXPLORATION OF 
REACTIVE SKELETON DYNAMICS 

Let us first detail how we provide initial conditions 
required to define Ti.sk and to solve Eqs. ([3"0|) . 



A. Initial conditions required to define reactive 
Skeleton dynamics 

It is obvious that the equations describing the 
late stages of binary black hole dynamics, given by 
Eqs. ([3H)l . require us to specify Hsk defined in terms of 
ipi(ip2,f,Pr, j) and •02 ("01 ; r,Pr,j)- This implies that we 
need to find ways to specify values for f,(f),p r ,j,ipi and 
ip2 at an initial instant, say t = 0. The initial conditions 
for f and <p are arbitrary and specify the initial radial 
separation and its associated angular position. 

The 2PN accurate expression for p r in terms of f in the 
ADM coordinates, displayed below, is our initial condi- 



tion for p r 

64 r\ /460 r] r) 2 \ 256 rjir 

f 1568 r? 2 158224 r, \ 

V 45 r 5 2835 r 5 / ' 1 ' 

The derivation of Eq. (|31|) requires following inputs and 
is done in three steps. First, we compute df/dt associ- 
ated with the quasi-circular inspiral to 2PN order with 
the help of Ref. 5] . Employing the 2PN accurate relation 
connecting j and f, relevant for PN accurate circular or- 
bits, we iterate the 2PN accurate Hamiltonian Equation 
for symbolically df/dt, written in terms of f,p r and j, 
and compute p r as a PN accurate polynomial in df/dt. 
In the final step, we invoke df/dt, accurate to 2PN radi- 
ation reaction order (and was computed in the first step) 
and obtain PN accurate expression for p r in terms of f, 
given by Eq. (|3"Tj) . 

The initial value for j is computed by solving dp r /dt = 
0, where dp r /dt is given by Eq. (|30c|) . It should be ob- 
vious that we need to specify r,p r ,ipi and ip2 to obtain 
numerically the value for j by imposing the right hand 
side of Eq. (|30c[) to zero. We provide initial guesses for 
ipi and tp2 by invoking the ingredients used to construct 
Brill-Lindquist potential function, associated with the 
Brill-Lindquist initial value solution for two black holes 
at rest [llj. This leads to the following expressions for 
ipi and "02 in terms of f and r\: 

^ =2-^f— + 2 V + -f + Af-*--^ (32a) 

^ = Lfi4Jl4 + l-l. (32b) 
1 2 4r 2 V r 4r 2 r 2 

We would like to emphasize that it is only during the eval- 
uation of an initial guess for j that we invoke Eqs. (|32p . 
Further, while performing various checks, we realized 
that the initial value of j is rather insensitive to the choice 
of p r . In other words, initial value for j does not change 
much even when we use p r — on the right hand side of 
Eq. (|3tB . 

We are now in a position to obtain initial values for ipi 
and ip2 that are required to solve Eqs. (|3"D|) . The initial 
values for ipi and ^2, associated with the above men- 
tioned f,p r and j values are numerically evaluated using 
the coupled algebraic equations for ipi and ip2, given in 
Eqs. (HI. 

From here onwards, we drop hat symbol appearing on 
our dimensionless quantities like TL^,f,p r ,j,Cj,t for the 
ease of presentation. 

B. Various aspects of the reactive Skeleton binary 
black hole dynamics 

Before we probe numerically the various aspects of 
the SAPE to model merging binary black holes, having 
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i] = 1/4, let us state that the way of solving Eqs. (130)) is 
quite different from the way one solves Hamiltonian equa- 
tions when the Hamiltonian is explicitly given in terms 
of r, p r and j. While dealing with the reactive 2.5PN 
accurate Hamiltonian system, given by Eqs. (|26[) . or the 
EOB scheme, one needs to invoke only once a numeri- 
cal scheme, like the 4 th order Runge-Kutta, to solve the 
associated couple differential equations, after specifying 
certain initial conditions for r, p r and j. However, a sim- 
ilar approach is not applicable for the reactive Skeleton 
dynamics, given by Eqs. ([50)) . This is because of the fact 
that the Skeleton Hamiltonian, given by Eq. ([2"2"| . implic- 
itly depends on r,p r and j and its explicit dependence 
on tpi and ip2 is interconnected. Therefore, the numer- 
ical solution to Eqs. (|30|) is determined in the following 
way. We have already detailed how to prescribe initial 
values for r, cf>,p r , j, ipi and ip% in the previous subsection. 
With the aid of these initial guesses, we invoke the 4 th 
order Runge-Kutta to solve Eqs. ([30)) and obtain values 
for r,p r ,(f> and j at the next time step St. These values 
are employed to evaluate the values of ipi and ip 2 at the 
instant St, with the help of the coupled algebraic equa- 
tions for ipi and ip2, given in Eqs. ([2"2"| . To obtain the 
values of the dynamical variables at the next instant, 25t, 
we use as initial guesses, the values of r,p r , 4>,j,ipi an d 
ip2 at the earlier step, and repeat the above detailed pro- 
cedure all over again till we reach the terminal point of 
our evolution. This is how we obtain numerical solution 
to the reactive Skeleton dynamics for binary black hole 
evolution, given by Eqs. (|30|) . In other words, we need 
to call the 4 th order Runge-Kutta at every time step t 
to obtain the values of r, p r , <f>, j, ip± and "02 at t + St. 
However, while solving the 2.5PN accurate Hamiltonian 
equations, Eqs. (|26[) or the EOB differential equations, 
one usually invokes a numerical scheme like the 4 th order 
Runge-Kutta only once. 

We are now in a position to extract (numerically) var- 
ious aspects of the reactive Skeleton dynamics for bi- 
nary black hole evolution. In this paper, we terminate 
the SAPE when u)(t) reaches its maximum value, w mx 
and let the initial value for r to be 10. Plots depict- 
ing temporal evolutions of various dynamical variables, 
appearing in Eqs. (|30|) . f° r a binary black hole having 
7] = 0.25 is presented in Fig. [TJ The parametric plot 
of x — r cos 4> versus y = r sin <f> clearly indicates that 
the dynamical evolution, predicted by Eqs. ([30)) , occurs 
along inspiralling quasi-circular orbits and it is rather 
difficult to pin-point the position of the last stable or- 
bit, which can clearly be defined using the conservative 
Skeleton Hamiltonian as done in Ref. [8]. Further, we ob- 
serve that both the binding energy, given by the Skeleton 
Hamiltonian, and j, the orbital angular momentum vary 
smoothly till we terminate the reactive Skeleton evolu- 
tion at oj — uJxax ~ 0.0896. The temporal behavior of p r 
clearly indicates that even in the neighborhood of u mx , 
the orbital motion is not substantially different from a 
quasi-circular inspiral. The region where p r increases 
sharply may considered to be the plunge phase. In Fig. [21 



we explore r(t),oj(t) and f 2 {t) under reactive Skeleton dy- 
namics. Plots of to versus t and cj versus r indicate faster 
evolution for lu after the late inspiral stage. The rapid 
decrease in r during the plunge phase is also evident in 
the r versus t plot. Interestingly, r 2 — (dr/dt) 2 remains 
small even during the plunge phase. However, it changes 
from ~ 10~ 3 to ~ 10~ 2 when t changes from ~ 588 to 
^ 627, which clearly indicate deviations from the adia- 
batic evolution along circular orbits of decreasing radii, 
expected during the quasi-circular inspiral. Another in- 
teresting aspect of the SAPE is that the reactive evolu- 
tion becomes slightly eccentric if we choose p r = at the 
initial instant t = and this is also displayed in Fig. [2] 
We would like to point out that residual eccentricities 
appearing in the puncture evolution in full General Rel- 
ativity were suppressed in a similar manner in Ref. [191 ]. 

The nature of the BBH evolution during the plunge 
phase, predicted by Eqs. ([50)) is probed in Fig. [3J First, 
we plot the oj 2 r 3 , the 'Kepler combination', against r. 
This plot is motivated by the observation in Ref. Q 
that for the reactive EOB evolution ui 2 r 3 , which remains 
close to unity during the quasi-circular inspiral, decreases 
rapidly during the plunge. We note that, similar to what 
is observed in the Numerical Relativity based evolution 
of punctures, the 'Kepler combination' varies even during 
the inspiral. Therefore, it is reasonable to state that the 
constancy of 'Kepler combination', observed in Ref. 0], 
is purely a consequence of the fact that the EOB scheme 
employs Schwarzschild-like coordinates. In Fig. [31 we also 
plot v 2 = w 2 / 3 versus t and Vq = r 2 ui 2 + f 2 versus t and 
observe that remains larger than vq throughout the 
Skeleton inspired inspiral. The above two plots indicate 
that in reactive Skeleton dynamics, it is not very attrac- 
tive to employ v = r u> or v = Vq in the expression for 
the far-zone angular momentum flux while incorporating 
radiation reaction on the conservative Skeleton dynam- 
ics. A plot of the 'quasi-circularity' condition dE = ujdj, 
valid for binary inspirals along a sequence of circular or- 
bits is also available in Fig. [3J The quantities dE and 
dj represent differences in the values of Hsk and j at t 
and t + St, while Q is the average value of w, given by 
Eq. (|30b[) , at t and t + St. We observe that the degree 
with which the above relation is violated keeps getting 
bigger as we approach a-> mx . However, the difference be- 
tween dE and Q dj is small even in the neighborhood of 
c<j mx indicating that the plunge motion is fairly close to 
'quasi-circular'. Therefore, in our opinion, the plots in 
Fig. [3] suggest that it is not very unrealistic to employ 
to estimate the effects of radiation reaction in our SAPE. 

Let us perform a first-order comparison between the 
GW phase evolutions based on SAPE and Numerical Rel- 
ativity (NR) for an equal mass binary during its late in- 
spiral without using any data from NR simulations. This 
is possible due to the observation, reported in Ref. [20| . 
that GW phase evolution associated with certain PN 
approximant, namely Taylor T4 at 3.5PN order, closely 
agrees with equal mass binary black hole NR simulations 
that last around 15 orbits prior to the merger [the accu- 
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mulated phase difference is iess than 0.05 radians over the 
30-cycle waveform] . The orbital phase evolution for equal 
mass binary black holes under TaylorT4 approximant at 
3.5PN order is obtained by numerically integrating the 
following two differential equations [201 ] : 



d(t){t) 

dt 
dx(t) 

dt 



,3/2 



(33a) 
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where 7 is the Euler gamma and x — uj 2 ^ 3 . In what fol- 
lows, we compare the orbital phase evolution prescribed 
by SAPE for 77 = 0.25 with <j)(t) originating from Eq. ([53]). 

There are two slightly different ways to compare </>sk(i) 
against 4>ta (t) after choosing a specific interval for x, de- 
fined by some minimum and maximum values, say x m i n 
and lEmax, for x. In the first approach, one lines up 
0Sk(i)j ^T4(i) and their time derivatives at x = x max and 
with the help of a shifted 0T4' (t) , given by Eq. (49) in 
Ref. [2(J, computes <ftsk(t) — 0T4'(*) when x = x m i n [ 
this is pursued in Refs. [2 if]. In the second prescription, 
one lines up (f>sk{t), 0x4 (i) an d their time derivatives at 
% — aWn a nd computes 0sk(£) ~ 0X4 (^) when x reaches 
a^max and this is what we list below and the differences 
between the two prescriptions in our comparisons are al- 
ways negligible. In Fig. 21 we plot 0sk(i) — 0X4 (i) in radi- 
ans and x(t) under SAPE and T4 approximant at 3.5PN 
order for various x intervals. For the top panel, the or- 
bital frequency range is from u — 0.025 to u> = 0, 05 and 
there are ~ 8.5 cycles. In the middle panel, we vary u> 
from 0.025 to 0.045 and there are about 8 cycles. Finally, 
for the bottom panel there are about 11.5 cycles and the 
uj interval is from 0.0222 to 0.045. Further, the fractional 
differences in are ~ 7%, 5.8% and ~ 6.1% in the above 
three cases and it increase sharply as we approach the 
plunge. 

It should be emphasized that we haven't employed 
any arbitrary parameters in our SAPE. Therefore, it is 
conceivable that by introducing flexible parameters, we 
should be able to obtain a modified reactive Skeletonian 
dynamics that prevents the above observed dephasing. 
This is inspired by what is advised in the modified EOB 
approach [4| and it will be reported in a separate publi- 
cation. 



IV. ROUGH PRESCRIPTION TO OBTAIN BBH 
COALESCENCE WAVEFORMS VIA SAPE 

We aim to provide, in this section, a rough estimate 
for the GW polarizations, having Newtonian accurate 



amplitude, associated with the BBH coalescence in our 
SAPE. This is achieved by terminating the SAPE when 
u>, given by Eq. (|30b[) . reaches its maximum value w mx 
an then matching there the relevant time derivatives of 
the radiative multipole moments associated with the late 
plunge phase to the corresponding 'ring-down' multipole 
moments, constructed using appropriate quasi-normal 
mode(QNM) contributions. Our construction of h™(t) 
and h+ (t) , i. e. the temporally evolving GW polarizations 
having Newtonian accurate amplitude, for the ring-down 
phase after the termination of the SAPE follows very 
closely what is detailed Sec. IV of Ref. [7[ . 

The GW polarizations are defined in the following way 



hs, = hi 



XX 

04> ' 



2\ oe <j><t> 



(34a) 
(34b) 



where and frJJ are the independent components 

of hj^- representing the radiation field emitted by iso- 
lated systems in asymptomatically flat space-time. The 
various angular components of /i^ T are defined using a 
spherical coordinate system, having orthonormal bases 
(f',6,4>), centered on the center-of-mass of the binary 
and whose azimuthal axis is aligned with the orbital 
angular momentum vector I. Following Ref. {22}, it is 
possible to obtain hj^ in terms of certain time deriva- 
tives of the radiative mass and current moments I lrn 
and S lm , irreducibly defined with respect to I, where 
m = — I, —1+1, .., +/ and T® 2 ' lm and T? 2 ' m : the so-called 
pure-spin tensor-spherical harmonics of electric and mag- 
netic types. To compute h^. , , we only require dominant 
contributions to hJ T given by 



TT 



G 



n = ^7 E Vl 2m (t-r'/c) T 

m=-2 



E2,2m 



(35) 



It turns out that at the Newtonian order, X 20 = I 21 = 
and contributions arising from m = — 2 are evaluated 
using the usual relation I 2 ~ 2 = I 2+2 , where * stands 
for complex conjugation. Therefore, only the following 
expressions are required to compute , 



Wi 22 = ~-Vw^v 



mc lo ' e 



(36a) 



T 



E2,22 
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1/2 



o 2i V> 



(1 + costf 2 ) 0® 



V128 7T, 

-0<g)0^ +2i cos #10 <g (/>+(/> 6\. (361)1 
With these inputs, we obtain 



h y — — if] cos i 



G 



" M W 2 / 3 sin2 0, (37a) 



^ = -2,(1 + 00^) (|£ 



w 2/3 cos2 0, (37b) 
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the orbital inclination angle similar to Eqs. ((37)) read 



where we identified i9 = 
and set ip = 0. 

Following Ref . Q , we aim to match as smooth as possi- 
ble (2)l piungc § iven in Ec l s - ©' at t = i m , t m being the 
time when the SAPE u(t) reaches its maximum value, 
to the corresponding ring-down ^I^ing consisting of cer- 
tain sum of decaying QNM modes. For simplicity, we 
restrict our attention to ^2^| constructed using the 
first two complex conjugated fundamental QNM modes. 
For this purpose, we only need to know the following 
conjugate pair of complex QNM frequencies, extractable 
from Ref. 12311 . 



(Gm\ 
— \ riuj{t) 8/3 cos 2 0(f) , 



(42a) 



3(V> 4 ) = -16 cosi ( ^ ) i)w(i) 8 / 3 sin 2 0(0 . (42b) 
\c z r' J 

The resulting plots for + (t) and the associated $t(ip4) 
and 9(^4) are displayed in Fig. [5] The plots clearly 
demonstrate a smooth matching to the QNM ringing 
phase in the SAPE. 



<t± = 0.08896 ± 0.37367 i. 



(38) 



Under these restrictions, ( 2 );r 2±2 during the ring-down 
phase is given by 



(2)2-2+2 _ Q-\ 



"ring ^0 

(for r = t-t m >0) 



(39) 



The two arbitrary complex coefficients Cq(X 22 ) and 
Cq (X 22 ) can be chosen so as to ensure not only that 
(2) ^pi 2 ungc(i = *m) agrees with ^l 2 l g {t = t m ) = 
Cq(2 22 ) + Cq(1 22 ), but also that the numerically com- 
puted time derivative of (2) I^ ungc (t), i.e. {3) I 2 ? ungc (t) 
agrees, when t = t m with ^J 22 ng = — a£ Cq(I 22 ) — 
o^Cq{1 22 ). This yields 



C+(Z 22 ) 



Co (I 22 ) 



- (3) T 22 ( f \ 
°20 x plungcA'v 



, (4) T 22 /.% 
' -^-plungcV'/ 

-1 

°20 ~ a 20 ' X 



+ (4) I 2 } (t) 
' -''plunge \ L / 



(40a) 



+ (3) T 22 ,.s 
°20 -'"plunge 



(40b) 



We would like to point out that while computing 
^-^piungeWi we include the dw/dt contribution, with the 
help of Eqs. (|A2[) . Using the above detailed prescriptions 
for ^I^ing ' ^ ^ s now fairly straightforward to obtain h™ + 
associated with BBH coalescence under SAPE that also 
includes contributions from QNM ringing phase. 

In Numerical Relativity, it is customary to invoke the 
Weyl scalar 1^4 to represent the gravitational wave con- 
tent of the dynamical space-time. In a suitable null 
tetrad and in the far-zone, one finds 

ipi = h + — i h x . (41) 

The 3t(ip4,) and 3(V , 4) having 'Newtonian' amplitude, 



V. DISCUSSION AND FUTURE DIRECTIONS 

This paper provides a semi-analytic prescription to 
model the BBH coalescence where we employ the punc- 
tures, employed by several Numerical Relativity groups, 
to model the non-spinning black holes. Our approach re- 
quires the conservative Skeleton Hamiltonian, derived in 
Ref. [H, representing two orbiting punctures in a wave- 
less truncation to the Einstein field equations. We in- 
clude the effects of radiation reaction in a transparent 
Hamiltonian framework that fully justifies the physical 
system we are trying to model. This is how we provide 
a prescription to implement our semi-analytic puncture 
evolution, termed SAPE, to model the merging binary 
black-holes. We also provided GW polarizations, real and 
imaginary parts of the Weyl scalar, having Newtonian 
amplitude, associated with the entire BBH coalescence 
in our SAPE. This is achieved by matching smoothly 
(2 ^piunge(*)' when w = w mx, to ( 2 >2?± 2 (t) , constructed 
with the dominant black hole quasi- normal modes. We 
observe that Eqs. ([50)) defining our semi-analytic punc- 
ture evolution display features and predictions qualita- 
tively similar to what are observed by various Numerical 
Relativity groups employing punctures to model black- 
holes. 

Efforts are being planned to compare the predictions 
of a suitably modified SAPE with that arising from the 
modified EOB approach. We have also initiated a pro- 
gram to include the spin effects into Hsk- Further, due 
to the fact that Hsk is exact while describing the mo- 
tion of a test black hole in the Schwarzschild spacetime, 
a refined version of the SAPE should be of some inter- 
est to GW physicists, excited by the prospect of LISA 
observing Extreme/Intermediate Mass Ratio Inspirals. 
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APPENDIX A: FORMULAE REQUIRED TO 
IMPLEMENT RIGHT HAND SIDES OF EQS. (J30|) 



A close inspection of Eq. (|22p reveals that 7igk depends implicitly on r, p r and j. Therefore, evaluating right hand 
sides of Eqs. ([30]) in the SAPE is rather involved. Using the theorem of implicit differentiation, dHsk/dp r may be 
computed using the following relation 

d Pr r ,j fixed I d Pr J \ dp r J J / I \d^2 J V^l/J 

and a similar relation holds for We follow a similar scheme to compute the right hand sides of equations that 

define d<p/dt and dp r /dt. 

During our matching to the QNMs, we require to compute ^-pi unge {t) arid this demands computation of dui/dt 
under the SAPE. From Eq. (|30bp . it is clear that u> requires us to compute dipi/dt and dfa/dt and we use following 
relations: 

#1 

dt 
(hp2_ 

dt 

«l 

We would like to emphasize Eqs. (|A2[) are only required during our matching to the QNMs. 



dr ) dt \ dp r J dt \ dj J dt 
dr J dt I dp r J dt I dj J dt 



(A2c) 
(A2d) 
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FIG. 1: Plots depicting various facets of SAPE, defined by Eqs. (|30[) . for a 77 = 0.25 BBH having initial orbital separation 
r = 10 [recall that we have dropped hat symbol appearing in the dimensionless variables defined by Eq. (|21|) ]. In these plots, we 
terminate the SAPE when u reaches its maximum value, u; mx ~ 0.0896. The Panel (i) provides a parametric plot of x — r cos cj> 
versus y — r sin0, indicating that the SAPE occurs along quasi-circular orbits. The temporal evolutions for the dimensionless 
binding energy, E = TCsk, and the orbital angular momentum j axe plotted in the panels (ii) and (iii). The plot for p r (t) 
suggests that even during the dynamical plunge, i.e. in the neighborhood of w mx , the orbital motion is not that different from 
a quasi-circular inspiral. 
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FIG. 2: Plots providing Lo(t), ui(r(t)), r(t) and r(t) 2 for BBH evolution under SAPE [the binary configuration is same as in 
Fig. [T] and all quantities are dimensionless] . The rapid changes in uj,r and r 2 occurring during the dynamical plunge is clearly 
visible and we terminate the SAPE when uj reaches its maximum value u; mx . Though r 2 = (dr/dt) 2 remains small near aj mx , 
its magnitude changes roughly fivefold during the plunge under SAPE. The black-hole binary evolution that employs Eq. (|3ip 
for the initial value of p r is displayed with thick line, while BBH evolution having p r = at the initial instant is displayed with 
dashed lines. The spurious eccentricity in SAPE is suppressed by using at t = the expression for p r , given by Eq. pip , and 
this is also reported in Ref. [TjJ. 
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FIG. 3: Plots displaying reactive evolutions of few interesting dynamical quantities and relations in the SAPE [the binary 
configuration is same as in Fig. [T]and all quantities are dimensionless] . The top panel plots v 2 , — ui 2 ^ 3 , vq = r 2 + r 2 uJ 2 . The 
strong r dependence of Vq makes it less attractive to characterize the orbital velocity compared to v u - In the middle panel, 
we plot the 'Kepler combination' (KC) uj 2 r s against r and a gradual decrease is clearly visible throughout the inspiral. The 
quasi-circularity (QC) condition (d E — lj d j) is displayed in the bottom panel and its smallness justify the use of in the 3.5 
PN accurate expression for dj/dt. 



15 




- — SAPE 

14 


/ 

/ 

J 1 


■ 1 , 






200 400 600 




800 



200 400 600 800 1000 1200 



0.075 
t 



200 400 600 800 1000 1200 



FIG. 4: Differences in the orbital phase evolutions associated with the SAPE and TaylorT4 approximant at 3.5PN order along 
with the associated x = lu 2 ^ 3 evolutions [recall that we have dropped hat symbol appearing in the dimensionless variables 
defined by Eq. (|21[) ]. From the top to bottom panels, the orbital frequency ranges are [0.025 — 0.05], [0.025 — 0.045] and 
[0.022 — 0.045] respectively. The fractional differences in (j> are ~ 7%, 5.8% and ~ 6.1% in the above three cases. A sharp 
increase in 0sk(i) — 4>Ti(t) as we approach the plunge is also observed. 




FIG. 5: Plots showing temporal evolution of scaled h x ,+ , 5ft(?/>4) and Q(ip4,) for equal mass binary black-hole coalescence in 
SAPE. We scale out Gm/c 2 r' from Eqs. (J3TJ) and (J42J) and let i = tt/3. The matching to the QNMs is based on Ref. 0] (see 
Section llVl and recall that t is identical to t, defined by Eq. (I2lj0 . 



